function [XX,V,VUnique,nVU] = interaction_expansion(X)


%Initialize
[n,p]=size(X);
V=bi2de(X);
VUnique = unique(V);
nT = numel(VUnique);
nVU = zeros(nT,1);
for ii = 1:nT
    nVU(ii) = sum(V == VUnique(ii));
end

%  "or" interaction expansion
W=zeros(n,2^p);
for j = 0:(2^p-1);
    Z = any(X(:,de2bi(j,p)==1)==1,2);
    W(:,j+1) = Z;
end
   
%  "and" interaction expansion
WW=zeros(n,2^p);
for j = 0:(2^p-1);
    Z = all(X(:,de2bi(j,p)==1)==1,2);
    W(:,j+1) = Z;
end

%Hadamard-Walsh Basis
WWW=zeros(n,2^(p-1));
for k=1:max(V);
    num1 = sum(de2bi(k-1,p));
    if num1 < 6 && num1 > 1;
        WWW(:,k) = (-1).^sum(X & (ones(n,1)*de2bi(k-1,p)),2 );
    end
end


% Wavelet Transform  (Not sure I'm doing this in the standard way)
%H = ConstructHaarWaveletTransformationMatrix(2^p);
%WWWW = H((V+1),:);

%Combine Bases
XX=WWW;






